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ABSTRACT 

Simulation design is the choice of locations in parameter space at which simulations are to be run and is 
the first step in building an emulator capable of quickly providing estimates of simulation results for arbitrary 
locations in the parameter space. We introduce an alteration to the "OALHS" design used by |Heitmann et &\\ 
(2006) that reduces the number of simulation runs required to achieve a fixed accuracy in our case study by a 
factor of two. We also compare interpolation procedures for emulators and find that interpolation via Gaussian 
Process models and via the much-easier-to-implement polynomial interpolation have comparable accuracy. A 
very simple emulation-building procedure consisting of a design sampled from the parameter prior distribution, 
combined with interpolation via polynomials also performs well. Although our primary motivation is efficient 
emulators of non-linear cosmological A/-body simulations, in an Appendix we describe an emulator for the 
cosmic microwave background temperature power spectrum publicly available as computer code. 
Subject headings: cosmic background radiation - cosmological parameters - methods: numerical - methods: 
statistical 



1. INTRODUCTION 

Many data analysis problems in cosmology require the re- 
peated computation of the power spectrum, or other sum- 
mary statistic, expected for a given point in the model pa- 
rameter space. These range from the very expensive, such as 
the calculation of weak lensing shear power spectra via ray 
tracing through A/-body simulations, to the relatively cheap, 
such as the calculation of cosmic microwave background 
(CMB) power spectra from Einstein-Boltzmann (EB) equa- 
tion solvers. Decreasing the computer resources and time re- 
quired for performing these calculations has benefits across 
this range of problems. 

The benefits of speed are perhaps most obvious at the more 
expensive end of the cost spectrum. As |Heitmann et al.] ( |2009] > 
have recently argued, a "brute-force" analysis to produce con- 
straints on cosmological parameters from "Stage IV" weak 
lensing shear power spectra 1 would take a 2048-processor 
machine 20 years. Techniques that can dramatically reduce 
the computer resource demands are not merely beneficial, but 
absolutely necessary. The "Cosmic Calibration" (CC) frame- 
work introduced by Heitmann et al. (2006) a nd developed in 
a series of papers (|Habib et al.||2007| |Heitmann et al. |2008 
Schneider et al.|20081| Heitmann et al.|2009[|Lawrence et al.| 



2009) has largely been motivated by this problem 



Even at the relatively cheap end there are benefits to in- 
creased speed. Although EB solvers are relatively fast (Seljak 
& Zaldarriaga|19 96 ), and can take as little as tens of seconds 



on a single processor, much effort has gone into making power 
spectrum calculation even faster. Early work exploited the 
angular-diameter distance degeneracy at small angular scales 
to reduce pre-compute requir ements (jTegmark & Zaldarriaga 
|2000| [Kaplingh atet al.|2002) . Later work used more sophis- 
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1 Such as from the Large Synoptic Survey Telescope (LSST) or a future 
space-based dark energy mission 



ticated training and interpolation procedures, and a lot of pre 
computing, to achieve very high accuracy at very fast post 
comp ute speeds ( |Fendt & Wandelt|2007bTa] |Auld et al.|2007 
2008 ). 

Following the terminology of Heit mann et al.| ( |2006| l, by 
emulator we mean a machine that can estimate the results 
of a power spectrum calculation, without actually doing the 
calculation. Typically emulator construction includes a pre- 
computing stage in which the calculation of the power spec- 
trum is performed at some finite set of locations in the param- 
eter space. The choice of this set of points (or the "simulation 
design") is a critical step in the emulator's construction. It 
can have a significant impact on the number of calculations 
required to achieve a given accuracy. The focus of this paper 
is on the development of simulation designs that minimize the 
computer resource demands of the pre-compute stage. 

The new design technique we introduce includes an appli- 
cation of the highly efficient "Orthogonal Array Latin Hyper- 
cube Sampling (OALHS)" which is part of the CC framework. 
The essential new elements are the choice of the parameter ba- 
sis in which to perform the OALHS and the volume reduction 
in the parameter space achieved by constructing designs in a 
hypersphere rather than a hypercube. This choice is informed 
by the (appropriately-weighted) use of information about the 
rates of change of the power spectrum with respect to param- 
eter variations. 

We also consider the interpolation procedures used for esti- 
mating the power spectrum at arbitrary locations in the param- 
eter space, from those pre-computed at the design locations. 
We compare results from the Gaussian process (GP) models 
used in the CC framework with those from fitting a second- 
order polynomial. In all cases we study we find the (much- 
simpler-to-implement) second-order polynomial fitting per- 
forms similarly to the GP interpolation. Their relative merits 
will vary from problem to problem however, and we speculate 
as to when GP will have advantages. We find that simulated 
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annealing is an enormously useful step to include for calibrat- 
ing the GP parameters. 

From our case study we draw conclusions that will be use- 
ful for application to other problems. Most importantly we 
demonstrate that the simulation design, rather than the choice 
of an interpolation method, is most critical for emulator con- 
struction. We also consider the choice of basis functions for 
reducing the dimensionality of the statistics to be interpolated, 
calibration of GP parameters for interpolation (and how this 
is affected by the choice of design algorithm), and the conver- 
gence of interpolation errors as the number of design points is 
increased. 

As a testing and development ground for our work on emu- 
lators we have chosen to work with CMB temperature power 
spectra. Although we are also motivated by weak lensing ap- 
plications, the relative speed of a Boltzmann-code calculation 
compared to an A^-body simulation makes the CMB tempera- 
ture power spectrum quite convenient for testing and develop- 
ment. Accurate calculation of non-linear matter power spectra 
is sufficiently expensive as to make their use for development 
of emulation strategies a practical impossibility. Heitmann 
|et al.| ( |2009| l circumvented this problem by developing and 
testing their design strategy using an approximate (but fast) 
method for calculating the matter power spectra (HALOFIT 
ISmith et al.|2003] ). 

While existing CMB emulators are very fast and sufficiently 
accurate, they do have very heavy pre-compute requirements. 
The PICO code relies on tens of thousands of evaluations. 
The heavy pre-compute needs make it difficult to extend an 
existing emulator's capabilities to include additional physical 
effects, such as isocurvature modes or completely new phys- 
ical models. Improved designs can greatly reduce these diffi- 
culties. As a demonstration we describe and validate a pub- 
licly available emulator in Appendix [A] called Emu CMB, 
that can calculate CMB temperature power spectra out to 
I = 5000 at sub-percent accuracy as a function of six cosmo- 
logical parameters. The development of the emulator, utiliz- 
ing the efficient design techniques we discuss here, only re- 
quired running an EB solver (CAMB Lewis et al. 20001 at 
100 design points. We also show in Appendix |A| the emulator 
performance for the CMB polarization power spectra, which 
still have large errors at low I. The large number of training 
points required for previous CMB emulators is in part driven 
by improving the accuracy for the low-^ polarization spectra 
(as well as covering a larger parameter space). 

The structure of this paper is as follows. We describe each 
aspect of our emulator in Section [2] We lay out the physical 
model for CMB power spectra that we use as a case study for 
the emulator construction in Section |2~T1 We consider 3 sim- 
ulation design methods that are described in Section |2T2| The 
basis decomposition of the simulation design runs and the 2 
interpolation methods we compare are given in Section 2.3 
and 2.4 We present the results of our case study in Section |3 
and draw conclusions about the performance and practical- 
ity of the different design and interpolation methods in Sec- 
tion |4] We present an example of a CMB temperature power 
spectrum emulator that meets the requirements for analysis of 
modern CMB experiments in Appendix [A] Code to construct 
emulators similar to the one in this section can be downloaded 
from |http : //www . emucmb . inf o| In Appendix [B] we 
compare the performance of quadratic polynomial interpola- 
tion in the matter power spectrum design of IHeitmann et al. 



TABLE 1 

Fiducial cosmological parameter values and \-a 

MARGINAL ERRORS DERIVED FROM THE FISHER MATRIX. 



Parameter 


Fiducial value 


WMAP5 error 


Planck eiTor 




0.959932 


0.0139280 


0.0048435 


too e* 


1.039648 


0.0033038 


0.0003797 


ln(10 10 A) 


3.049634 


0.0475892 


0.0346602 




0.1079211 


0.0069907 


0.0017480 


Ub 


0.022490 


0.0006217 


0.0001885 


T 


0.086429 


0.0175253 


0.0162973 



2. EMULATORS DESCRIBED 

There are three steps to constructing an emulator. 

1. Choose a "simulation design," i.e., the set of points 
in the parameter space at which the summary statistic 
(henceforth assumed to be a power spectrum) will be 
calculated. 

2. Perform a reduction in the number of dimensions of the 
power spectrum, via an (incomplete) mode decomposi- 
tion. 

3. Specify an interpolation procedure to allow one to es- 
timate the power spectrum that one would have calcu- 
lated at any point in the parameter space. 

Here we review how these three steps a r e implemented in th e 
"CC" work of IHeitmann et al.| pOOoTl; |Habib et al.| pOUT]); 
Heitm ann et al.| ( |2008) >; [Heitmann et al.| (|2009|); |Law rence 
|et al.| ( |2009| ) and also with PICO ( |Fendt & Wandelt|2007b|a| l 

2.1. CMB Temperature Power Spectrum For Our Case Study 

The very low expense of CMB temperature power spectrum 
calculations makes them convenient to use for testing and de- 
velopment of new emulation techniques. For this reason, we 
pursue emulation of a CMB temperature power spectrum cal- 
culator for our case study. 

We set ourselves the goal of building an emulator with suffi- 
cient accuracy for analyzing the temperature power spectrum 
constraints that are expected from Planck 2 , and to do so in 
the simplest possible manner with the fewest possible "sim- 
ulations." To inform the emulator construction, we assume 
prior constraints on cosmological parameters of similar qual- 
ity to those from WMAP5 ( |Nolta et ak||2009| ). We consider 



the following six cosmological parameters: 

0= (n s ,lOO0*,ln(lO lo A) ,u> c ,u> 6 ,t) 



(1) 



where n s is the scalar spectral index, 9* is the angular size 
of the sound horizon at last scattering (in radians), A is the 
amplitude of the primordial power spectrum defined at the 
pivot scale k = 0.002 Mpc" 1 , lo c = Vl c h 2 is the dark matter 
density, cut, = flbh 2 is the baryon density, and r is the optical 
depth. Our fiducial values for these parameters are given in 
the first column of Table [T] and are derived from a Markov 
Chain Monte Carlo (MCMC) chain using the WMAP5 likeli- 
hood 3 . We use a Fisher matrix to compute constraints on 
these six cosmological parameters given a signal power spec- 
trum computed with CAMB and a model for the noise power 
spectrum, 



(2009) with the GP interpolation errors presented in their pa- 
per. 
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FIG. 1 . — Error models for the CMB temperature power spectrum used in 
the Fisher matrix predictions for WMAP5 (black) and Planck (yellow). 

with noise weights w no i se = 1.43 x 10 14 for WMAP5, w n0 ise = 
4.7 x 10 16 for Planck, beam smearing cr beam = FWHM/2.355, 
and FWHM = 0.00378 for WMAP5 (at 90 GHz) and 0.0021 
for Planck (at 143 GHz). 

We further include in our Fisher matrix a prior on a combi- 
nation of optical depth and primordial power spectrum ampli- 
tude similar to what is provided by WMAP5 polarization data. 
Specifically, assuming the quantity Ae~ 2r has negligible error 
(since this combination is determined very well from the tem- 
perature power spectrum), we impose a Gaussian prior with 
width: 



a (At 2 ) cr(e 2T T 2 ) 2(t + 1)ct(t) 



At 2 



(3) 



with ct(t) = 0.017 from WMAP5 4 , which comes from the 
combination of polarization data with temperature data and 
is the main contribution of the polarization data to parameter 
constraints. In Figure [T] we show the size of the power spec- 
trum error bars for our WMAP5-like observations, and those 
for our Planck-like observations. In Figure [2] we show the 
corresponding two-dimensional Fisher matrix contours. 
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0.0219 0.0225 




FIG. 2. — 2-D marginal parameter constraints predicted by the Fisher 
matrices for Planck (red, solid) and WMAP5 (black, dashed) noise and beam 
models. 



We were sufficiently satisfied with the accuracy and sim- 
plicity of the emulator we developed, that we are releasing a 
version of it as a publicly available code, called Emu CMB, 
described in Appendix |A| There is another fast and accurate 
code available (PICO IFendt & Wandeitl|2007b|a| l that has 
some advantages. Most notably, PICO is a more complete 
emulator that provides predictions for polarization as well as 
temperature power spectra and includes models with non-zero 
spatial curvature. Emu CMB's advantages are (1) its signif- 
icantly simpler algorithm that has made it very easy for us to 
make available implementations in multiple languages with 
no need for external libraries, (2) the simple training algo- 
rithm can be reimplemented quickly (using our public code if 
desired) for revised CMB models, and (3) it extends to higher 
maximum multipole moment (£ max = 5000). The latter qual- 
ity is especially important for analysis of high-resolution data, 
such as from the Atacama Cosmology Telescope or the South 
Pole Telescope (SPT). 

The design runs for this paper were all calculated using 
CAMB (2008 June release [Lewis et al.|2000| > with accuracy 
settings accuracy_boost = 4, l_accuracy_boost 
= 4, 

l_sample_boost = 8, and without including gravita- 
tional lensing 5 . With this accuracy setting, numerical inte- 
gration errors are negligibly small. We modified CAMB to 
extract the power spectra only at multipoles where they are ac- 
tually computed, skipping the default interpolation to all mul- 
tipoles. For a maximum multipole of I = 3000, this gives 448 
power spectrum values using the default multipole spacing in 
CAMB. 



4 |http : //lambda ■ gsf c ■ nasa . gov/product/map/dr3/| 

|params/lcdm_sz_lens_wmap5 ■ cfm| 



2.2. Simulation design construction 

Here we review two design procedures: that used for CC 
and that used for PICO. We also introduce our new design 
procedure. We compare results from these three design pro- 
cedures in Section [3] 

The CC design approach is "Orthogonal Array Latin Hyper- 
cube Samphng^OT^ |Morris & M itchell 
[1995] |Leary et al.|2003| |Carnell|2009f 7lt begins with defin- 
ing the boundaries of a hypercube in a scaled parameter space 
- all parameter dimensions rescaled to fit in an interval be- 
tween and 1. A grid in this hypercube is then defined. In 
a two-dimensional parameter space the «</ design points are 
assigned to the grid in such a way that each design point is 
the only occupant of its row and the only occupant of its col- 
umn. In a 3-dimensional space, the occupation of site {i,j,k} 
means that sites {l,j,k} are empty for all I i, sites {i,l,k} 
are empty for all / j and sites {i,j,l} are empty for all 
I y k. The generalization to n-dimensional spaces is straight- 
forward. The final step is to then set the exact location of each 
design point within its grid cell. This is done using various 
criteria that generally aim to maximize the distance between 
points in any lower-dimensional projection. The fact that the 
points are well spread out, in many different ways of quantify- 
ing "spread out", improves interpolation to non-design points. 

The design for PICO is constructed from samples from a 
prior probability distribution for the parameters and hence we 

5 The effects of gravitational lensing are included in EMU CMB, since it 
is an important effect to include for analysis of real high-resolution data. 
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TABLE 2 
Parameter range for OALHS 
design 



Parameter 


Minimum 


Maximum 


n s 


0.9042 


1.0156 


1000* 


1.0264 


1.0529 


ln(lO'°A) 


2.8592 


3.2340 


w f 


0.07996 


0.1359 


Ub 


0.02000 


0.02498 


T 


0.01633 


0.15653 



call this design procedure PS for "Prior Sampling". In partic- 
ular, the PICO design is a draw of samples from a converged 
MCMC run. The design is motivated by the fact that the sam- 
ples are drawn from exactly the region of the parameter space 
where we desire accurate emulation (i.e., the region with non- 
negligible prior probability). 

The parameter bounds used in our OALHS design are given 
in Table [2] The bounds have a width of eight times the 
marginal Fisher matrix errors for the WMAP5 noise model, 
centered on the fiducial point given in Column 2 of Table [T| 
Given these bounds and the number of design points, rid, we 
draw a realization of an OALHS using the maximinLHS 
function of the R Ihs package ( |Carnell |2009[ l, which maxi- 
mizes the minimum distance between the design points within 
the Latin hypercube. All two-dimensional projections of 
a 100-point OALHS design are shown in the upper triangle of 
panels in Figure [3] 



0.02 0.08 0.14 



.4ft, 



ln(10 10 A) 



0.020 0.023 

FIG. 3. — Location of 100 simulation design runs in the 6-dimensional 
parameter space. Upper triangle is for the OALHS design and the lower 
triangle is for the OALHSFS design (black triangles) and the PS design (red 
crosses). 

Both the OALHS and PS design procedures have their ad- 
vantages. The OALHS design fills the parameter space in an 
efficient way, providing well-spaced anchor points for inter- 
polation. The PS design, on the other hand, concentrates de- 
sign points only in the region of parameter space where we are 
actually interested in evaluating the emulator. We now present 
a design procedure that captures the advantages of both. 

If pe is the number of input parameters to our simulation, 
we can compute a Fisher information matrix (based on the 
data error distribution determining the parameter prior) from 
the summary statistic of the simulations we wish to emulate 



using 2pe + l additional simulation runs; that is via pe two- 
point numerical derivatives of the summary statistic and one 
fiducial point used to calculate the summary statistic covari- 
ance model. We use this Fisher matrix to rotate the input pa- 
rameter axes to a less-correlated set of input parameters. We 
then build an OALHS design in this de-correlated parameter 
space, but reject all samples of the OALHS design that are 
outside of a hypersphere of constant probability (at quadratic 
order) centered on the fiducial parameter location. We call 
this design OALHSFS for OALHS samples inside a "Fisher 
sphere". We also use this decorrelated space for performing 
the interpolation of the simulation design runs. 

The two-dimensional projections of a 100-point OALHSFS 
design constructed from a Fisher matrix using the WMAP5 
noise model are shown as the black triangles in the lower 
panels of Figure [3] The design points are located within the 
hypersphere of radius 4 in the decorrelated parameter space 
(note that each parameter has unit variance in this space). We 
chose a radius of 4-cr because this encloses over 98% of the 
volume of a six-dimensional Gaussian distribution, which we 
judged to be sufficient for our emulator. (Note that this is cov- 
ering 98% of the volume of the prior distribution. The frac- 
tion of the posterior distribution covered will be much larger 
.) The points of a 100-point PS design selected as random 
samples from an MCMC chain using the CMB power spec- 
trum with the same WMAP5 noise model are overplotted as 
red crosses for comparison. It is readily apparent that the 
OALHSFS design points are confined to the region of high 
prior probability as defined by the MCMC samples. 

By choosing design points inside a hypersphere in a de- 
correlated parameter space, we can achieve large reductions 
in the volume of the space that must be covered by the design 
points as the dimensionality of the parameter space increases. 
In three dimensions, the volume of a sphere embedded in a 
cube, such that the sphere's diameter is equal to the length 
of a side of the cube, is nearly half the volume of the cube. 
In six dimensions, the hypersphere has a volume nearly 12 
times smaller than that of the surrounding hypercube. In prin- 
ciple then, if the interpolation error in our emulator is largely 
determined by the mean distance between design points, we 
should be able to achieve similar interpolation precision using 
12 times fewer design points in an OALHSFS design as in an 
OALHS design in a six dimensional parameter space. Or, for 
a fixed number of design points, the OALHSFS design should 
have smaller interpolation errors than the OALHS design. 

As shown in Figure [3] the PS design gives similar volume 
savings over the OALHS design as does the OALHSFS. How- 
ever, the PS designs do not consistently cover the design space 
and can lead to very large interpolation errors in regions where 
the design points are sparse 

The Fisher matrix approximates the shape of the likelihood 
about its peak at quadratic order in the cosmological parame- 
ters. Our OALHSFS design will fail to encompass a region of 
constant probability if higher order corrections to the shape of 
the likelihood are significant. The PS design does not suffer 
from this problem. In practice, it will be possible to measure 
the accuracy of the OALHSFS probability contours by com- 
paring with a PS design as we have done in the lower panels 
of Figure [3] If necessary, the OALHSFS design can be ad- 
justed to allow for different ranges along the different axes in 
the "de-correlated" parameter space. That is some axes might 
extend for, e.g., ±4er while some might extend for ±12cr, etc. 
(where a is the marginal error determined from the Fisher ma- 
trix). 
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In addition to the prior parameter ranges, the parameter 
ranges along each axis of the simulation design could also 
be informed by the curvature of the mode amplitude surfaces 
to be interpolated. In the case that a design axis has a range 
much shorter than the response surface curvature scale along 
a given coordinate axis the design will provide an estimate of 
the curvature scale. This could result in, e.g., the response sur- 
face being modeled as less smooth than it really is. Informa- 
tion about the curvature of the response surfaces is contained 
in the second derivatives of the power spectra with respect to 
the cosmological parameters (while the Fisher matrix uses the 
first derivatives of the power spectra). The range of the design 
axes could then be increased in directions where the curvature 
is very small (e.g., using the mode amplitude with the small- 
est curvature to set the axis length). On the other hand, if the 
response surface is known to have an exceptionally high cur- 
vature along a given axis, the density of design points could 
be increased to accurately sample the high curvature. Because 
the emulators for our case study and Emu CMB work to an 
acceptable precision we have not explored these ideas further. 

It is always possible that when inferring parameter con- 
straints from a data set via MCMC, the MCMC chain will 
move some steps out of the region of parameter space covered 
by the emulator design. If the MCMC chain spends very little 
time in the region outside the design, it may possible to fall 
back on evaluating the actual Boltzmann code for analyzing 
the CMB. For applications with more expensive simulations 
(such as A^-body simulations) the emulator's can predictions 
can be extrapolated to get an initial MCMC chain that can 
then be used to inform the construction of a new design. 

2.3. Dimensional Reduction 

To reduce the dimensionality of the function to be interpo- 
lated, both the CC and PICO emulators perform a principle 
component decomposition of the n y x n^ matrix of simulation 
design runs y\ = l(£+l)/(2ir)C' g where i = I, rid an d ^ takes n y 
values in the range (2,£ max ). The principal components that 
contribute least to the variance over the design are discarded. 
We have found that keeping the 20 most-important modes is 
more than sufficient for our purposes, i.e., the errors made 
from neglecting the other min (n Y , nj) modes are insignificant. 
We quantify the errors from the truncated mode decomposi- 
tion in two ways. First, we have verified that using 20 PC 
modes allows us to reconstruct the power spectra at all design 
points in several OALHS designs to within 0.1%. Second, we 
used the Fisher matrix to estimate the cosmological parameter 
biases that would be induced from the sy stematic errors in the 
reconstructed power spectra (see Huterer & Takada 2005 for 
details on the Fisher matrix methods). For = 20, the distri- 
bution of parameter biases had first and third quartiles < 10% 
for all six of our cosmological parameters. This means that 
in the case of perfect interpolation of the power spectra, the 
inferred parameter constraints using our emulator would have 
biases less than 10% of the 1 — a marginal errors for both the 
OALHS and OALHSFS design s. 

Following IHeitmann et al | ( p006| ; |Habib et ai] ( |2007l ); 
Schneider et al. ( 2008) > we first subtract the (^-dependent) 
mean of the design runs from each design power spectrum 
and then scale the result by a single number so that the com- 
bined entries of our n y x nj matrix of design runs have vari- 
ance one. Denoting this centered and scaled matrix r/, we then 
perform a singular value decomposition, rj = UBV 7 " where U 
has dimension n y x p (p = min(n y ,nd)) with U r U = I p , V has 



dimension nj x p with V r V = I p , VV T = I nj , and B (p x p) 
is a diagonal matrix of singular values. Finally we define the 
matrix of basis vectors, $ = — !=UB and w = \ fruV J so that 

iw r w = L. 

Retaining only the p^ < p most significant columns of <& 
we can write the scaled power spectrum as a sum over p„ 
modes 

(4) 



,x=l 

where is the entry of the matrix $ in row £ and column /i, 
and w^(6) is the ^th (parameter-dependent) mode amplitude. 

In Figure|4]we show the first four PC basis functions for the 
power spectra in the three different design types, PS, OALHS, 
and OALHSFS. Because of the very different volumes of pa- 
rameter space covered, application of the PC algorithm to 
the design runs yields very different basis functions for the 
OALHS design versus the PS and OALHSFS designs. This 
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FIG. 4. — The basis functions corresponding to the first four most signif- 
icant principal components in the decomposition of our three types of 100- 
point designs. 

is a crucial difference because the basis functions must ac- 
curately describe the design spectra at all multipoles in order 
to meet the accuracy targets over the large dynamic range we 
are simulating. Note that modes 1 and 2 derived from the 
PS and OALHSFS designs have several oscillations, but are 
much smoother when derived from the OALHS design. This 
difference can be seen later in the interpolation errors from 
the different designs (Figure [7]) where the OALHS design has 
larger, and more oscillatory, errors at large multipoles indicat- 
ing that the missing oscillatory structure in the OALHS design 
basis functions leads to missing structure in the reconstructed 
power spectra. In Appendix |A| we address the choice of basis 
functions again when we extend the emulator for the CMB 
temperature power spectrum to I = 5000. 

2.4. Interpolation 

We compare two methods for interpolating the mode am- 
plitudes Wfj,{ff) between design points in the cosmological pa- 
rameter space. 

For the first interpolation method we again follow the CC 
framework and model w^(6) as independent GPs for each 
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value of ijl. We re fer the reader to |Habib et aL] ( |2007| ); |Heit-| 
|mann et aT7| ( |2QQ9| > for details about the construction of the GP 
model for the mode amplitudes. In summary, for a fixed set 
of points Of (i = 1 , . . . , «</), the amplitudes w^ithetcii) are mod- 
eled to have a multivariate Gaussian distribution with mean 
zero and a parametrized covariance structure that determines 
the smoothness of the interpolation. The covariance param- 
eters are calibrated from a likelihood model conditioned on 
the simulation design runs. Because the GP framework infers 
distributions for the covariance parameters rather than best fit 
values, the error in the GP fit and interpolation can be propa- 
gated through to the final cosmological parameter constraints. 

Our second interpolation method is to fit a /^-dimensional 
linear or quadratic polynomial to the design runs. Compared 
to the GP model, this has the advantage of being very fast to 
compute, but the disadvantage of not propagating any inter- 
polation error. For small numbers of design points, the or- 
der of the polynomial is limited (in order to obtain a well- 
conditioned fit) and this can further limit the accuracy of the 
polynomial interpolation. We therefore expect GP interpola- 
tion to outperform global polynomial fits when the curvature 
of the mode amplitude surfaces (in any parameter direction) 
exceeds the maximum polynomial order that can be fit with 
nj design points. A caveat is that GP models with our chosen 
covariance will also yield unacceptably large interpolation er- 
rors for extremely rapidly varying surfaces unless the number 
of design points is significantly increased. So there is a gen- 
eral limit on the smoothness of the mode amplitudes as func- 
tions of cosmological parameters for such high-dimensional 
parameter spaces and limits on the practical number of design 
runs. 

Figure [5] shows the first and fourth mode amplitudes in the 
PC decomposition of the power spectra in the OALHS 100- 
point design as a function of n s with all other cosmological 
parameters fixed at their central values. The range of n s covers 
the marginal 4-a error as derived from the Fisher matrix with 
the WMAP5 noise model. The first mode amplitude in this 
case is a very smooth function of n s , while the less significant 
modes have a higher degree of curvature. (Note that the size 
of the residuals in the mode amplitudes does not obviously 
translate into power spectrum errors because of the centering 
and scaling of the power spectra done before performing the 
mode decomposition.) 

For our case study we have found the step of calibrating the 
GP parameters to have intensive computatational and human- 
labor requirements. With the GP covariance parametrization 
used in the CC framework there are p^ variance parame- 
ters and x po correlation parameters. To calibrate the GP 
model we must then find the maximum of the likelihood in 
this high-dimensional parameter space. We used MCMC to 
find the vicinity of the GP likelihood peak, but found the most 
difficult step to be initializing the chains in regions of non- 
negligible probability. A simulated annealing algorithm with 
a constant cooling schedule turned out to be quite effective in 
getting our GP calibration chains started. Because the calibra- 
tion of the GP model can be a large computational obstacle in 
constructing an emulator we highly recommend the use of a 
simulated annealing or similar algorithm. We do not, how- 
ever, know of any algorithm that can replace the human labor 
of assessing the convergence of the multiple MCMC chains 
required to optimize the GP model. We demonstrate the per- 
formance of our simulated annealing algorithm for calibrating 
the GP interpolation parameters in Figure [6] 

We also found the choice of the proposal distribution for 



the GP correlation parameters in the MCMC algorithm to be 
crucial to efficient calibration of the GP model. We used in- 
dependent uniform distributions for each correlation param- 
eter, but initially allowed the width of the uniform proposal 
to be adjusted according to the variance in the previous 1000 
steps in the chain. Because this breaks the detailed balance 
requirement for convergence of the MCMC chain, we fix the 
widths of the proposal distribution after a preset number of 
chain steps (of order 10 5 ). 

Because the evaluation of the GP likelihood at each chain 
step requires inversion of an «</ x matrix, the GP interpo- 
lation method becomes unpractical as becomes very large. 
This is unfortunate because the interpolation accuracy should 
improve as n.4 increases! We therefore expect GP models to 
provide a preferable interpolation method for small and 
polynomial interpolation to be preferred for large n^. For our 
case study, we actually find that polynomial interpolation pro- 
vides smaller interpolation errors for a wide range of n^. We 
demonstrate the performance of polyno mial interpolation for 
a different case study matching that of Heitmann et al. ( 2009) > 
in Appendix [B] 

Our C++ code for calibrating GP models for interpola- 
tion (built with the Scythe Statistical library |Pemstein et al. 
|2007 1 is available for download at |http : / / www . emucmbT 
linfol 

3. EMULATOR RESULTS COMPARED 

3.1. Designs Compared 

In Figure [JJ we show the fractional interpolation error for 
100 test spectra using the three design methods under consid- 
eration. The test spectra were computed at points in our six- 
dimensional cosmological parameter space drawn from the 
posterior distribution for our WMAP5 noise model. Each in- 
terpolation method was trained on 100 design points. (Note 
that we always show the expected value of the emulator; we 
do not compute the random draws that are part of a complete 
GP prediction.) At these test points, the OALHSFS design 
has interpolation errors approximately two times smaller than 
the OALHS design at all multipoles and for both interpolation 
methods. Using quadratic polynomial interpolation, the PS 
design has comparable interpolation errors to the OALHSFS 
design, but has very poor performance using GP interpola- 
tion. As we mentioned in Section |Z2j the OALHSFS design 
is built in a parameter-space volume that is ~ 1/12 the size 
of the OALHS volume. If the increased density of simulation 
points is the dominant difference between the OALHSFS and 
OALHS designs, then Figure [7J gives some indication of how 
the interpolation accuracy improves with the increased den- 
sity. However the interpretation is not straightforward as the 
shape of the designs (spherical versus cubic) may also be im- 
portant. That is, the design points outside the Fisher sphere 
in the OALHS design could, in principle, be helpful in reduc- 
ing the interpolation accuracy, but we have not pursued the 
separation of these effects further. The main result that the 
OALHSFS design gives improved performance is sufficient 
for our purposes. We show similar results, but for = 50 in 
Figure [8] The interpolation errors using quadratic polynomi- 
als have similar relative performance between the three design 
types as in Figure [7J The GP interpolation results, however, 
now show smaller errors using the PS design than using the 
OALHS design. The OALHSFS design using GP interpola- 
tion generally gives smaller errors than the OALHS design, 
but there are a few test points where the OALHSFS design 
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FIG. 5. — The first and fourth mode amplitudes for the OALHS 100-point design as a function of n s , with all other comological parameters fixed at their central 
values. The range of variation in n s is given by the marginal 4-cr errors as inferred by the WMAP5 Fisher matrix (i.e., the total width of the OALHS design along 
the n s axis). The gray lines show GP draws conditioned on the design runs using the maximum-likelihood GP parameters. The blue lines show the mean of the 
GP draws. The red lines show a polynomial fit that is quadratic order in the 6 cosmological parameters. 
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FIG. 6. — Trace plot of the likelihood in an MCMC chain using simulated 
annealing to calibrate the parameters of a GP model for design interpolation. 
This particular chain was calibrating a PS design. 

gives worse performance for several multipoles. These results 
are consistent with our expectation that at a fixed number of 
design points nj the reduction in the mean distance between 
design points in the OALHSFS and PS designs (see Figure [3} 
should improve the interpolation accuracy. Indeed we do see 
an improvement by a factor of ~ 2 in the interpolation accu- 
racy. However, the poor performance of the PS design with 
100 design points and using GP interpolation demonstrates 
that the PS design is a less reliable algorithm than OALHS- 
based designs for building an emulator using GP interpola- 
tion. This could be because the PS design leaves regions 
of parameter space very sparsely covered compared to the 
OALHS designs (especially in higher dimensions). 

Next we investigate how the accuracy gains of OALHSFS 
over OALHS at fixed nj translate into a reduction in nj at 
fixed interpolation accuracy. In Figure [9] we show the quar- 
tiles and maxima of the average fractional interpolation error 
in a band from / = 1500 to I = 2000 (that is, we averaged the 



absolute value of the error over I = 1500 to 2000) versus nj 
for 100 test points taken from a PS design. We plot the inter- 
polation errors for both the OALHS and OALHSFS designs 
using GP interpolation. Comparing the solid and dashed lines 
at fixed interpolation accuracy (i.e., along a horizontal axis) 
we see the OALHSFS design requires ^25-70% fewer de- 
sign points than the OALHS design (with larger savings at 
low rid). One possible explanation for this trend is that some 
of the design points in the OALHS design are not significantly 
contributing to the reconstruction of the mode amplitude sur- 
faces in the (high-probability) region where the test points are 
actually located. When nj is small, every design point should 
be important for reducing the interpolation errors, so wasting 
even one point in the OALHS design would degrade the em- 
ulator. This demonstrates a main point of the paper that the 
design space volume reduction offered by the OALHSFS de- 
sign is significant in minimizing rid to achieve a fixed error 
tolerance. 

In order to understand when the errors in the truncated 
mode decomposition of the power spectra are significant rel- 
evant to the interpolation errors, we plot the fractional er- 
rors in the reconstructed power spectra as functions of 
at each design point (i.e. without any interpolation) in Fig- 
ure [10] The reconstruction errors are rapidly decreasing func- 
tions of the number of retained PC modes, as expected. For 
the choice used throughout this paper of p^ = 20 for rid > 50, 
the errors from the truncated mode decomposition are ade- 
quate for our error tolerance of 10~ 3 for both the OALHS 
and OALHSFS designs. However, from Figure [9] we can see 
that with pp. = 20 the mode reconstruction errors are non- 
negligible for the rid = 100,200 designs, which would ex- 
plain why the total fractional power spectrum errors (mode 
reconstruction + interpolation) do not steadily decrease for 
n d > 100. 

3.2. Interpolation Methods Compared 

In Figure [TTJ we compare the mean interpolation errors in 
the band 1500 < I < 2000 using the GP model and a quadratic 
polynomial fit to the OALHSFS and PS designs as functions 
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FIG. 7. — Fractional power spectrum interpolation errors using 3 different 100-point simulation designs: PS (left), OALHS (center), OALHSFS (right). The 
interpolation is performed using GP (top) and quadratic polynomial (bottom) fits to the 20 most significant principal components in each design, and is tested 
against 100 points in an additional PS design. The dark grey bands show the minimum and maximum errors out of the 100 design points for each / while the red 
(inner) band shows the errors within the first and third quartiles. The central black line shows the median interpolation error. 



0.004 




500 1000 1500 2000 2500 3000 



500 1000 1500 2000 2500 3000 



500 1000 1500 2000 2500 3000 



FIG. 8. — Same as Figure|7]except for nj = 50. 

of rid. In a ll of the designs considered polynomial interpola- tion provides smaller interpolation errors than GP interpola- 
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FIG. 9. — Fractional error on the interpolated power spectra (using GP 
interpolation) averaged over 1500 < £ < 2000 as a function of the number 
of design points for 2 designs: OALHS (black - solid) and OALHSFS (red 
- dashed). Thick lines show the maximum of the absolute value of the two 
quartiles of the error distribution over the 100 test points. Thin lines show the 
maximum absolute error. 
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FIG. 10. — Fractional errors in the power spectra reconstructed from the 
truncated mode decomposition averaged over 1500 < £ < 2000 as a func- 
tion of the number of Principal Components retained. Our standard choice 
throughout the paper is = 20. The power spectrum errors are computed 
at each design point for each n c [ considered in Figure [9] (See the Figure [9] 
caption for the definitions of the quantile labels). Note the points have been 
offset horizontally for clarity. 

tion for small rid- For the OALHSFS design, the interpolation 
methods give comparable errors for large rid. The PS designs 
have sporadic error ranges, indicating the large variation in 
the covering of the design space test points with different PS 
design realizations. 

The (first or second order) polynomial interpolation per- 
forms so well because the response surfaces of each mode am- 
plitude are so smooth as functions of the six active cosmologi- 
cal parameters. The GP models also require smooth response 
surfaces, but we expect the low-order polynomial fits to be 
worse than the GP fits for a given correlation length in the 
response surface as the correlation length is decreased. This 
is because the covariance structure in the GP model has more 
flexibility than a (global) second order polynomial. There- 
fore, for some applications the fast-to-compute polynomial 
interpolation will no longer suffice. The relative performance 



of different interpolation methods will always be problem- 
specific. However, we explore designs for the matter power 
spectrum in Appendix [B] and show that the polynomial inter- 
polation is again comparable to the GP interpolation demon- 
strated in Heitmann et al. ( |2009| l for similar designs for most 
of the dynamic range in the simulated power spectra (except 
for the smallest scales). So, for both the CMB and matter 
power spectrum, an emulator constructed using an OALHSFS 
design and quadratic polynomial interpolation has similar per- 
formance to an emulator using an OALHS design and GP in- 
terpolation. As the number of design points increases, both 
interpolation methods should converge to arbitrary accuracy 
as higher-order polynomial fits become well-conditioned. 

4. CONCLUSIONS 

We have demonstrated an improved method for generat- 
ing simulation designs for interpolation based on selecting 
OALHS inside a spherical (rather than cubical) volume in a 
whitened cosmological parameter space defined by the Fisher 
matrix. Our improved design algorithm (OALHSFS) offers a 
significant improvement in design efficiency over the OALHS 
algorithm; giving a factor of ~ 2 reduction in the number of 
design points (rid) required to reach a fixed emulator accuracy. 
We also compared the OALHS-based designs with designs 
constructed by drawing samples from the prior distribution for 
the cosmological parameters (which we call prior sampling or 
PS designs). Because both PS and OALHSFS designs con- 
centrate design points in the volume of parameter space with 
significant prior probability, they both yield smaller interpo- 
lation errors for a fixed rid than OALHS designs when rid is 
large. When rid is small, the PS designs give poor coverage 
of the parameter space; giving erratic interpolation errors de- 
pending on the particular PS design realization. This is not 
surprising because it is exactly the problem that OALHS de- 
signs are meant to solve. However, we take the superior per- 
formance of both the PS and OALHSFS designs in the large 
rid limit (rid = 100 for our case study in six-dimensional pa- 
rameter space), as evidence that the OALHSFS design retains 
the advantages of the OALHS design while also capturing the 
benefit of limiting design points to regions of significant prior 
probability. 

We compared the GP interpolation methods used in pre- 
vious literature of the emulator framework for cosmological 
simulations with simple, global polynomial interpolation over 
the design space. We find that second order polynomial inter- 
polation gives interpolation errors that are only slightly worse 
than a GP model with optimized correlation parameters. This 
is a useful discovery because the GP models can be compu- 
tationally challenging to calibrate while a polynomial fit via 
least-squares minimization is very fast. For PS designs we 
found it is not always possible to fit a GP model, while sec- 
ond order polynomial interpolation still gives errors similar 
to those using an OALHSFS design. For applications where 
low-order polynomial interpolation does not provide suffi- 
cient precision, we note that the GP model calibration done 
via MCMC can be significantly hastened by using a simulated 
annealing algorithm to search for regions of high likelihood in 
the GP parameters. 

Making use of the simplicity of polynomial interpolation, 
we are releasing code to generate emulators of the CMB 
power spectra whose computational limitation is only the time 
it takes to run the Boltzmann code to generate spectra at the 
design points (available at http : / /www . emucmb . inf o\. 
We demonstrate one such practical emulator in Appendix [A] 
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FIG. 1 1. — Fractional error on the interpolated power spectra averaged over 1500 < t < 2000 as a function of the number of design points using different 
interpolation methods: quadratic polynomial (black - solid) and GP (red - dashed). The left panel shows the interpolation errors in the PS designs while the right 
panel shows the OALHSFS designs. Thick lines show the maximum of the absolute value of the two quartiles of the error distribution over the 100 test points. 
Thin lines show the maximum absolute error. 



for the TT power spectrum reaching £ > 5000 and including 
lensing contributions. This emulator uses 100 design points to 
achieve interpolation errors less than 0.6% in the power spec- 
trum for all I. Python and IDL scripts as well as a module for 
CosmoMC ( [Lewis & Bri dle 2002) for generating predicted 
power spectra from this emulator are available at the same 
URL. 

Because of the similar performance of the two interpola- 
tion methods we investigated, we conclude that the simula- 
tion design is the most significant component of the emulator 
framework introduced in Heitmann et al.| ( |2"006| ); |Habib et al.| 
( [2007] ); [HeTtmann et al/ 



2009). The mathematically consis- 



tent GP framework becomes important when it is necessary to 
propagate the errors in the design interpolation or when opti- 
mizing the interpolation for a given OALHS-base d design. 
As mentioned by other recent work ( |Heitmann et al.|2006 



Habib et al. 2007 ; Heitma nn et al.|2009] l the simulation emula 
tor technology becomes most interesting for emulating com- 
putationally expensive simulations such as those needed for 
predicting tracers of the matter power spectrum (e.g., weak 
lensing or galaxy spectra). In this case A^-body simulations 
of the dark matter distribution can take days or weeks to run, 
which puts severe constraints on the number of design points 
that can be considered. In addition to the matter power spec- 
trum, the need for predicting the outputs of suites of simu- 
lations over cosmological parameter space has recently been 
recognized for understanding the deviation from universal- 
ity of the halo mass function (Ma nera et al.||2009] l, the dis- 
tribution o f progenitor m asses in the conditional mass func- 
tion fNeistein et al.|2010|, computing the 1 -point probability 



distribution function of the Lya flux ( Viel et aL]|2009| l, and 
mapping t he input space of semi-an alytic models of galaxy 
formation (Bow er et al.|20 09 2010). We expect the methods 
developed here to be useful for the first three applications. 
However, the problem of mapping the input space for simula- 
tions with a large number of parameters is qualitatively differ- 
ent. Our design based on the Fisher matrix assumes that the 
prior probability distribution in the input parameter space is 



unimodal and that the maximum of the prior probability dis- 
tribution is known a priori. The success of the interpolation 
methods discussed in this paper also assume that the ranges 
of input parameters are sufficiently small that the outputs are 
(quite) smooth functions of the input parameters. Both of 
these assumptions are broken when mapping the input space 
for semi-analytic galaxy formation models so this should be 
considered a separate problem for emulator construction (see 
Bower et al. 2010| for how OALHS designs are useful for 
such problems). 

Recently, An gulo & White| ( |2010| l have shown that statis- 
tics such as the power spectrum (as well as mass functions 
and merger trees) can be obtained at different cosmological 
parameter values using appropriate rescalings of a single N- 
body simulation. If this method proves to be successful in fur- 
ther tests, it could alleviate the strict bounds on the number of 
design points that are feasible for statistics derived from cos- 
mological A^-body simulations. In this case, the advantages of 
GP interpolation methods would become less useful (namely 
the ability to optimize the GP model to achieve minimal in- 
terpolation errors and to propagate the interpolation errors). 
Alternatively, using the methods of | Angulo & White] ( |2010| l 
could make it feasible to construct emulators of more sophis- 
ticated statistics derived from light cones requiring many N- 
body simulations for a single point in parameter space. We 
plan to pursue this in future work. 
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APPENDIX 



EMU CMB: A CMB TT POWER SPECTRUM EMULATOR FOR CURRENT EXPERIMENTS 

In this section we present a CMB TT power spectrum emul ator suitable for u se in parameter estimation algorithms from 
Planck 6 , ACT 7 , SPT 8 , or similar experiments. We use CAMB ( |Lewis et al.|2000| to generate CMB TT power spectra with the 
same six cosmological parameters as in the main body of the paper*, but out to £ max = 5238 and including the lensing distortion of 
the power spectrum. We have not concerned ourselves with modeling uncertainties such as those due to recombination. Instead, 
we regard this emulator as proof that a practical tool for current data analyses can easily be constructed. Emulators for revised 
models can be quickly constructed using the same algorithm. We provide scripts for this purpose as well as details of the exact 
CAMB settings that generated the results presented in this section at |http : / /www . emucmb . inf o| 



Building the emulator 

We use an OALHSFS design with = 100 covering the WMAP5 4-cr error bounds as determined by our Fisher matrix 



described in Section 2.1 The design runs for Emu CMB use the 2009 October release of CAMB (different from the case 
study in the main body of the paper) with the "curved correlation function" lensing implementation including nonlinear cor- 
rections calculated by HALOFIT (Challino r & Lewis||2005| l. We again use accuracy settings accuracy_boost = 4 and 
l_accuracy_boost = 4, but build the emulator from the power spectra as output by CAMB at every multipole rather than 
extracting only the uninterpolated power spectra values. We set l_sample_boost = 4 to determine the number of multipole 
values where the power spectra are actually computed. See |http : / /www . emucmb . inf o] for the complete CAMB parameter 
files. 

Also in contrast to the main body of the paper, we perform the dimension reduction of the design runs on y\ = 
log (£(£+ 1)/ (2tt) C' e ) (i= 1 , rid). Performing a PC decomposition on the logarithm of the power spectra yields basis functions that 



Al 



better describe the lensing variations at I > 3000 across the design points as shown in Figure . 

As in the main text, we fit global second order polynomials to each of the first 20 mode amplitudes to perform the interpolation 
over the design space. Using only the first 20 PCs allows reconstruction of the power spectra at the design points to less than 
0.1% accuracy for £ < 100 and to less than 0.01% accuracy for £ > 100. 



Emu CMB accuracy 

We validated the emulator by running CAMB at an additional 100 points in a new realization of an OALHSFS design. (Note 
the validation plots in the main body of the paper use points from a PS design.) The interpolation errors for the spectra at these 



100 test points are shown in Figure A2 For all test points, the emulator is accurate to less than 0.6% out to £ = 5238. The errors 
are less than 0.2% for £ < 3000, which is consistent with Figure]?] 

Note that while Emu CMB uses significantly fewer training runs of CAMB than PICO, unlike PICO this Emu CMB example 
does not predict polarization power spectra. We have verified with simple tests that interpolation of the low-i? parts of the EE and 
TE CMB power spectra is more challenging than the TT spectra demonstrated here. In addition, many of the training runs for 
PICO were required to include cosmological models with nonzero spatial curvature, which we have not considered 10 . 



http : / /www . esa . int /esaMI /Planck/ index . html 



http : //www .physics . prince ton . edu/act 



http : / /pole . uchicago . edu| 



We revert to the more common definition of the primordial amplitude defined at the pivot scale k = 0.05/iMpc 
C. Fendt private communication. 
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FIG. Al. — Left: TT power spectra at each of the 100 design points for EMU CMB. The mean of the design runs is shown in red. Right: basis functions for 
the 5 most significant principal components of the 100 design runs log(£(£ + l)/(27r)Q). 



MATTER POWER SPECTRUM 

In this section we compare the p olynomial inter polation errors for OALHS designs for th e dark matter power spectru m to the 
GP interpolations demonstrated in Heitmann et al. ( 2009|). F ollowing Heitmann et al. (2009) we use the HALOFIT (Smith et al. 



2003) 1 algorithm as provided in CAMB ( Lewis et al.|2 000) to generate nonlinear dark matter power spectra as functions of the 
five cosmological parameters: il c h 2 , Qbh A , n s , w, a%. The reduced Hubble parameter h is determined from uj c and u)\, assuming 
spatial flatness and a fixed angular size of the sound horizon at last scattering (6* = 7r/302.4 steradia ns as defined in|Hei tmann 
|et al.|2 009). The prior parameter ranges defining the boundaries of the OALHS design are the same as|Heit mann et a l.|(|2009| r 

The fractional interpolation errors for three designs with «</ = 10,37, 100 are shown in Figure Bl Heitm ann et al.| ( |2009) > built 
an OALHS design with = 37 a nd d emonstrated sub-percent interpolation errors using GP interpolation over the same range in 
wavenumber as shown in Figure Bl Using quadratic polynomial interpolation, the errors are worse, but are still under 2% for 
all wavenumbers. Increasing «</ to 100 points reduces the interpolation errors to sub-percent levels for all but a few test points 
in the design. For the matter power spectrum over this design space, the GP interpolation developed by Heitm ann et al.| ( |2009) > 
therefore requires three times fewer design points to achieve the same accuracy as a naive quadratic polynomial fit. However, 
we note that the OALHSFS design demonstrated in Section 3.1 can achieve the same reduction in the required number of design 
points to reach a fixed accuracy while using the much simpler polynomial interpolation. 

We take these results as confirmation that the good performance of polynomial interpolation in the CMB designs in the main 
body of the paper is not limited to a single physical example. More generally, the success of a simulation emulator depends 
strongly on the simulation design (including the form of the mode decomposition of the power spectrum) while optimized GP 
interpolation can also significantly improve the performance of a given Latin-hypercube-based design (in contrast to a PS -based 
design). 
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FIG. B 1 . — Fractional error in the interpolated matter power spectra u sing polynomial interpolation for 3 OALHS designs with different tij. The parameters 
and widths of the OALHS designs match those in Heitmann et al. 1(2009) . For = 10 the mode amplitudes in the decomposition of the power spectra are fit with 
linear polynomials. For nj = 37, 100 the mode amplitudes are fit with quadratic polynomials. 



